Finite-size effects in transport data from Quantum Monte Carlo simulations 
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We have examined the behavior of the compressibihty, the dc-conductivity, the single-particle 
gap, and the Drude weight as probes of the density-driven metal-insulator transition in the Hub- 
bard model on a square lattice. These quantities have been obtained through determinantal quantum 
Monte Carlo simulations at finite temperatures on lattices up to f 6 x 16 sites. While the compress- 
ibility, the dc-conductivity, and the gap are known to suffer from 'closed-shell' effects due to the 
presence of artificial gaps in the spectrum (caused by the finiteness of the lattices), we have es- 
tablished that the former tracks the average sign of the fermionic determinant {{sign}), and that a 
shortcut often used to calculate the conductivity may neglect important corrections. Our systematic 
analyses also show that, by contrast, the Drude weight is not too sensitive to finite-size effects, being 
much more reliable as a probe to the insulating state. We have also investigated the influence of 
the discrete imaginary-time interval (Ar) on (sign), on the average density (p), and on the double 
occupancy (d): we have found that (sign) and p are more strongly dependent on Ar away from 
closed-shell configurations, but d follows the Ar^ dependence in both closed- and open-shell cases. 
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PACS numbers: 71.10.Fd 71.30.+h 71.27.+a 73.63.-b 

I. INTRODUCTION 

Metal-insulator transitions (MIT) are still a topic of 
intense activity.'^' In clean systems, an otherwise metallic 
system can become an insulator through the opening of 
a gap in the spectrum due to electronic repulsion; they 
become what are known as Mott insulators.^ Alterna- 
tively, band insulators correspond to systems in which 
the valence band is completely filled, even in the absence 
of repulsive interactions. When the on-site energies are 
different (but regularly distributed), due to, say differ- 
ent atomic species, electrons may become trapped: in 
this case the system is a charge-transfer insulator. In 
addition, in the presence of disorder the system may be- 
come an insulator as a result of electrons being unable 
to diffuse throughout the lattice; i.e., they may undergo 
an Anderson localization transition. One clear experi- 
mental signature of the insulating state is a vanishing 
conductivity as the temperature is decreased. However, 
from the theoretical point of view, and in the context of 
Quantum Monte Carlo (QMC) simulationa^nZ! in partic- 
ular, detecting an insulating state is not always straight- 
forward. First, one necessarily deals with systems of fi- 
nite size, hence with gaps in the spectrum which may be 
of the same magnitude as the ones responsible for the 
insulating behavior. These gaps occur at filling factors 
corresponding to 'closed shells', and give rise to atyp- 
ical behavior in several quantities of interest; further, 
these closed-shell effects, which are readily seen in the 
non-interacting case, can persist in the presence of in- 
teractions (see below). Secondly, QMC simulations are 
plagued by the 'minus-sign problem','' which precludes 
the study of several low-temperature properties of the 
system as the electronic density is varied continuously. 
And, finally, in spite of the wide variety of quantities at 
our disposal to probe a MIT, such as the compressibil- 
ity, the dc-conductivity, the single-particle gap, and the 



Drude weight ,'Si to name a few, they yield conflicting in- 
formation in some cases, the origin of which is still not 
fully understood. For instance, under somewhat restric- 
tive conditions^ti^ the dc-conductivity can be calculated 
in a convenient way, without resorting to analytic contin- 
uation of imaginary- time QMC data to real frequencies, 
which may be a delicate matter fSEll however, in the case 
of the Hubbard model, for some particular combinations 
of lattice size and electronic densities (away from half 
filling), the conductivity behaves as if the system were 
insulating, which casts doubts on whether the conditions 
are really met, or if it is a manifestation of 'closed-shell' 
effects, or both. In the case of homogeneous versions of 
well studied models, one may be able to generate data for 
many different lattice sizes for a given electronic density 
('minus-sign problem' permitting); in this way, a trend 
with system size can be established, and any deviation 
from it should be readily identified. However, this may 
not be the case of s ystem s with an overlying structure, 
such as a superlatticej ^^ l ^^ l a checkerboard lattice, or even 
in the presence of staggered on-site energies (the ionic 
Hubbard model) .E^HUl 

Our purpose here is to shed light into these discrep- 
ancies, and to compare different approaches to detect a 
MIT from QMC data; as a by-product, we will also estab- 
lish a connection between the behavior of the compress- 
ibility and the infamous sign problem of the fermionic 
determinant. The layout of the paper is as follows. In 
Sec. [IT] we introduce the Hubbard model, and outline the 
computational approach used. In Sec. |III| we discuss the 
predictions from the electronic compressibility, when the 
effects of closed shells manifest themselves as a major 
finite-size effect. Section IIVI is devoted to finite-size ef- 
fects on the dc-conductivity and the density of states, 
as obtained through an inverse Laplace transform of the 
current-current correlation function and the the single 
particle Green's function respectively; in this Section 



we also provide numerical estimates for the errors in- 
volved when the dc-conductivity is calculated setting the 
imaginary-time r = j3/2, where, as usual, /3 = 1/T, in 
units such that the Boltzmann constant is unity. In Sec. 
[V] wc discuss the Drude weight in detail, and show that 
it docs not suffer from closed-shell effects. The single- 
particle excitation gap is considered in Sec. |VI[ and we 
find that it suffers from the same closed-shell effects as 
the other probes of the insulating state, apart from the 
Drude weight. In Sec. |VII| a systematic study leads to a 
connection between the sign of the fermionic determinant 
and the compressibility; we also discuss the influence of 
the imaginary-time interval on some of the data. And, 
finally, Sec. |VIlil summarizes our findings. 



II. MODEL AND CALCULATIONAL DETAILS 

The simplest model to capture the physics of Mott 
insulators is the repulsive Hubbard model, which is char- 
acterized by the Hamiltonian 

1\ / 1 

2 



+UY^ \n,^^- j [n,y 



M 



E"- (1) 



where, in standard notation. 



is the fcrmion dcstruc- 



tion operator at site i with spin a =t,i, rii^ = cl^c^^, 
and rii — Ui-^ + Ui^. We only consider nearest-neighbor 
hopping (indicated by {i,j)) on a two-dimensional Lx L 
square lattice, and work in the grand-canonical ensem- 
ble; the chemical potential /z is tuned to yield the desired 
density p = '^i{ni)/N, where N = L'^ is the number of 
lattice sites. The hopping parameter t sets the energy 
scale, so we take t = 1; throughout this paper, we have 
considered the weak- to intermediate coupling regime, 
[/ < 4, for which size effects are more severe. 

We use d etermi nant quantum Monte Carlo (DQMC) 
simulationa^iEEM] ^q investigate the properties of the 
Hubbard model. In this method, the partition func- 
tion is expressed as a path integral by using the 
Suzuki- Trotter decomposition of exp(— /3?^), introducing 
the imaginary-time interval At. The interaction term 
is decoupled through a discrete Hubbard-Stratonovich 
transformation,^'^ which introduces an auxiliary Ising 
field. This allows one to eliminate the fermionic de- 
grees of freedom, and the summation over the auxiliary 
field (which depends on both the site and the imaginary 
time) is carried out stochastically. Initially this field is 
generated randomly, and a local flip is attempted, with 
the acceptance rate given by the Metropolis algorithm. 
The process of traversing the entire space-time lattice 
trying to change the auxiliary field variable constitutes 
one DQMC sweep. For most of the data presented here, 
we have used typically 1,000 warmup sweeps for equili- 
bration, followed by 4,000 measuring sweeps, when the 
error bars are estimated by the statistical fluctuations; 
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FIG. 1. (Color online) Electronic density vs. chemical po- 
tential: in (a) for the free case, at /3 = 16, and for different 
linear lattice sizes L; in (b) for the L = 10 lattice and differ- 
ent interactions U . The horizontal dashed lines highlight the 
specific density (p = 0.42) at which one plateau appears for 
the L = 10 lattice. 

when necessary, the data were estimated over an aver- 
age of simulations with different random seeds. Typi- 
cally, we have set At = 0.125, but often data were also 
collected for At = 0.0625, just to confirm that system- 
atic errors are indeed small; further, for some quantities 
we have also performed extrapolations towards At — > 
from up to eight distinct values of At. One should also 
keep in mind that since we do not use a checkerboard 
breakup of the lattice, our equal imaginary-time data for 
[/ = are exact, so that they do not depend on the 
imaginary-time discretization; the T-dependent quanti- 
ties result from sampling even for [/ = 0, but the statis- 
tical errors are negligible in this case. Wit h the updating 
being carried out on the Green's functions,EHIllat the end 
of each sweep we have at our disposal both equal-'time' 
and T-dependent quantities, which we discuss in turn. 



III. ELECTRONIC COMPRESSIBILITY 



Let us first consider the electronic compressibility, 
K = p~^^dp/d^. Being a direct measure of the charge 
gap, it may be used to detect insulating phases; a major 
computational advantage is that it is a local quantity, 
thus fluctuating very little within the DQMC approach. 
In Fig. fTla), the density p is plotted as a function of the 
chemical potential, for different lattice sizes, for the free 
case, [/ = 0, and at a fixed temperature. If taken at 
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FIG. 2. (Color online) Total energy iJ as a function of 
electronic density p for the L = 10 lattice at T = in the 
noninteracting limit ([/ = 0). 

face value, plateaus in the p x fJ- curves would be iden- 
tified with incompressible phases, and hence with insu- 
lating regions. However, a closer look reveals that both 
the width of the plateaus, as well as their positions, are 
strongly dependent on the finite system sizes used. Given 
that ior U — the system is certainly metallic for all 
densities, the presence of these plateaus can be traced 
back to gaps in the energy spectrum of the noninteract- 
ing Hubbard model on a finite square lattice, which is 
given in the usual way by E ^ J2q<qp{p);a^i^)' '^i*^ 
e(q) — —2t{cosqx + cosQy), where q_F{p) is the Fermi 
wave vector for the density p. In Fig. [2] the total en- 
ergy is shown as a function of the electronic density for 
a 10 X 10 lattice: the energy gaps do not have the same 
magnitude, and one should notice, in particular, the gap 
at p — 0.42, which is quite large in comparison with the 
ones between levels with E < —2. This gap appears as 
a plateau in the data for the 10 x 10 lattice in Fig. llTa), 
indicated by the horizontal dashed line. The existence of 
this 'gap' is a manifestation of what is referred to as 'the 
closed-shell problem' and is characteristic of the finite- 
ness of the lattice. It should be stressed that such effects 
are still present when the interaction is switched on, at 
least up to intermediate values of U: from Fig. [ijb), we 
see that the gap moves towards smaller values of ^ as C/ 
is increased, though without any noticeable decrease in 
magnitude; in what follows we illustrate further conse- 
quences of these closed-shell effects. As the lattice size is 
increased, the gaps become smaller, and the plateaus in 
the electronic density become narrower, until they com- 
pletely vanish in the bulk limit, L — > oo. For this reason, 
from now on we will refer to these plateaus as pseudo- 
insulating states. The use of the compressibility to locate 
insulating regions must therefore be supplemented with 
thorough analyses of the robustness and the width of the 
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FIG. 3. (Color online) DOS spectrum for different lattice 
sizes at p = 0.42, and for [/ = (a) and f/ = 2 (b). The 
Fermi energy {to = 0) is shown as a dashed line. The error 
bars represent statistical errors from different realizations. 



plateaus with system size and temperature. 

IV. CONDUCTIVITY AND DENSITY OF 

STATES 

The optical conductivity and the density of states 
(DOS) are other probes of the insulating state which 
are worth discussing in depth; this is especially in or- 
der, given that the use of the shortcut to calculate 
the dc-con ductiv ity (see below) has been increasingly 
widespread,l^2l2Sl even beyond QMC.^Sl 

First we recall that the simulations yield imaginary- 
time quantities, such as the real-space single-particle 
Green's function, 

G(r^i-j,T) = (ci^(r)c]jO)), < r < /3, (2) 

and the current-current correlation functions, 

A(q,T) = (j,(q,r)j,(-q,0)), (3) 

where jx (q, r) is the Fourier transform of the 
'time'-dependent current-density operator, jj.(i, r) = 
e^^J2;(i)e"^^, with 



jxii)=itJ2{4+st,a% 
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FIG. 4. (Color online) Optical conductivity from inverted 
Laplace transform (see text) at (a) half filling, p — 1, and (b) 
for p — 0.42, at given U and inverse temperature, for different 
lattice sizes and At. The error bars represent statistical errors 
from different realizations. 



Now, the fluctuation-dissipation theorem yields^Zl 



A(q = 0,T)- 
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ImA(q = 0,w), (5) 



and hnear response theory impHei 
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IniA(q = 0,a;) = ujRea{uj); 
similarly, we have^ZISIl 



G(r = 0,r) 
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The calculation of cr(a;) and N{uj) is then reduced to 
numerically invert these Laplace transforms at a given 
temperature. Here we employ an analytical continuation 
method,'^ through which the conductivity and the DOS 
can be obtained for the whole spectrum w.E^ While there 
has been some debate over which type of analytic contin- 
uation me thod is best suited to perform these Laplace 
transforms ,'^^'22] our purpose here is not to perform a 
systematic study of the outstanding issues; instead, we 
adopt one of the procedures'^ to extract estimates for 
a{uj) which, in turn, will be used to test the trends in the 
calculation of (Jdci &s discussed below. 

In Fig. [3] we compare the DOS at density p — 0.42 
for the free and interacting cases, obtained through the 
method described in Ref. [151 It is clear that irrespective 
of the value of U, the DOS vanishes at the Fermi energy 
for L = 10, while being non-zero for L = 6 and 12. Figure 
[4] shows the optical conductivity for the interacting case 



(U — 2), calculated with the same inversion method,^ 
both at half filling, and for p = 0.42. We see that, while 
at half filling the insulating behavior is apparent for all 
system sizes [adc{T) — limi^^o o-(aj, T) — > 0], for p = 0.42 
one would be led to identify an insulating behavior if 
only data for a 10 x 10 lattice were available. One should 
also note that data for both Ar = 0.125 and 0.0625 are 
the same, within error bars. The origin of this 'false 
insulating' behavior can therefore be traced back to the 
closed shell effects discussed above, though here cr(w) is 
particularly affected by the large gap required to add an 
electron to the 'closed shell' of 42 electrons; an analogous 
problem occurs at the closed shell density of p = 86/144 
for the 12 X 12 lattice (not shown). 

In addition to suffering from the closed shell problem, 
the inversion procedure adopted'^J can be very costly in 
computer time, due to the need of very small error bars 
in the data for A. An alternative methocP^ to obtain 
CTdc(T) consists of setting r = {3/2 = 1/2T in Eq. (|5|, and 
assuming cr(w) admits a Taylor expansion near w = 0; the 
integral can, in principle, be carried out term by term in 
the surviving even powers of w, and we get 



^dc(T)«a("j(r) 

plus higher order terms, with 
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Note that if one wants to take a^J (T) as an approxima- 

tion for o'dc(T'), cr^J (T) must be small; this should occur 
if the temperature is low enough, and the frequency de- 
pendence of the conductivity is smooth, i.e., ii T <^ fl, 
where fl sets a small energy scale of the problem.'^ While 
it is hard to assess a priori if this condition is satisfied, in 
the present case we have data for cr(w) at our disposal, for 
several temperatures; this allows us to calculate t^^ (T), 
and check the errors involved in neglecting it in Eq. (|8|. 

For p = 1, we see from Fig. Mja) that cr^J{T) rises as 
the temperature is lowered [(red) circles], but eventu- 
ally bends down at some temperature, consistently with 
a^^J — > as T — > 0; in a generic situation, in which QMC 
data for these lowest temperatures were not available, 
one could be misled to state that the system is metal- 
lie. However, when a^J (T) is included [(blue) triangles 
in Fig.pTa)], the conductivity acquires the correct steady 
decrease with decreasing T ^ 0.15; this shows that higher 
order terms may indeed be crucial at temperatures not 
so low. Figure p[b) shows that for p = 0.42, cr^J steady 

decreases as T decreases, which is suggestive of insulating 

(2) 
behavior; the inclusion of data for cr^J {T) does not revert 

this trend. Since for other densities the metallic behav- 
ior is unequivocal [see, e.g., Fig.ISJc)], one concludes that 
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FIG. 5. (Color online) Temperature dependence of the dc- 
conductivity [circles, zeroth order in uu, and triangles, up to 
second order; see Eqs. (|8|-(10l], and of the Drude weight 
(squares; see Sec. Wh, for [/ = 2, on a 10 x 10 lattice, and 
for different electronic densities. The error bars for adc are 
due to the averaging process, while those for D/e^ are due to 
extrapolations towards aj,„ — > 0. 



the spurious effect for p = 0.42 is yet another manifesta- 
tion of the closed-shell density. We have found that these 
overall features are also present for U ~ A; in particular, 

(2) 

the contribution of cr^J {T) when p — 0.42, though signif- 
icant, is again not sufficient to yield a metallic behavior, 
thus confirming that the false insulating state is indeed 
a closed shell effect. 

From this analysis we conclude that extreme care must 
be taken when examining limT->.o f^^ (T) to indicate 
whether the ground state is metallic or insulating; in ad- 
dition, while the overall trend may be captured (away 
from closed-shell densities), attempts to fit experimental 



data with cr^^ (T) should lead to error, if the tempera- 
tures involved are not too low. 



V. THE DRUDE WEIGHT 



We now discuss the Drude weight, D, defined through 



lim Re (j{uj,T) = D6{uj) + a,^g{uj), (11) 
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FIG. 6. (Color online) Current-current correlation function 
A(q = 0, uJm) at half-filling p — 1.0 , as a function of ujm/'2nT, 
where Um is the Matsubara frequency, at a fixed inverse tem- 
perature /3 — 16. The solid symbols denote {—kx). In (a) 
the on-site repulsion is kept fixed and the data correspond 
to different linear lattice sizes; in (b) data are for a 12 x 12 
lattice, but for different values of U. The error bars represent 
statistical errors from different realizations. 



where crrog(i^) is the regular (or incoherent) response. 
Approxim ants to D are readily available from QMC sim- 
ulations a^^E2l 

^^^^ = [(-fc,)-A(q = 0,i..„)], (12) 

where a;„i = 2rmTT is the Matsubara frequency, and {kx) 
is the average kinetic energy of the electrons per lattice 
dimension. The Drude weight is then given by 



hm D„, = Hm D(T) = D. 

r,m-i-0 T^O 



(13) 



In actual calculations, both limits should be taken 
through extrapolations of sequences of low-temperature 
frequency-dependent data Dm,(T)^ finite-size effects 
and finite-Ar effects must also be taken into consider- 
ation when analyzing the data. 

FigurepTa) illustrates how the uniform current-current 
correlation function at half filling depends on the Mat- 
subara frequency, with both /3 and U fixed, for different 
system sizes. While {—kx) is hardly dependent on the 
system size (see solid symbols in Fig. l6|, the same does 
not hold for A{0,Ldjn)- Nonetheless, approximants to the 
Drude weight, as given by Eq. (12 1, do indeed approach 



zero with growing linear lattice size L, as it should for an 
insulating state. Figure [6jb) displays the same quantity, 
now for a fixed system size, but for different values of U; 
we see that as m — ?► 0, £),„ — >■ for [/ 7^ 0, while Dm 
approaches a non-zero value for U = 0. 
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FIG. 7. (Color online) Same as Fig. |6) but for p = 0.42; in 
(b) data are for a 10 x 10 lattice. 
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FIG. 8. (Color online) Drude weight approximants as a func- 
tion of temperature for p = 1.0 for a L = 12 lattice for differ- 
ent values of U. The error bars result from uncertainties in 
the extrapolations ujm/2TTT — ^ 0; see text. 



Data for p — 
We see that D,, 



0.42 and 

7^e2 [Eq. 



U 



(121 



2 are shown in Fig. [7] 
for the 10 x 10 lattice 



does not show any false insulating behavior, as it did for 
other quantities: in (a) the difference between {—k^) and 
lim(^^_j.o A(0,a;,„) does not display a significant change 
with lattice size, while in (b) the data show that the 
closed shell problem does not manifest itself over a wide 
range of values of U. 

In order to extract more quantitative data, we adopt 
the following procedure: For fixed L, C/, and /3, we plot 
D„^ as a function of m = uJm/i-KT, and extrapolate to 
771 — >■ with the aid of a parabolic fit to the data for 



FIG. 9. (Color online) Size dependence of the (normalized) 
Drude weight at a fixed, finite temperature, /? = 16, for differ- 
ent densities: p = 0.42 (circles), p — 0.66 (squares), and p = 1 
(triangles). Empty, half-filled, and filled symbols respectively 
correspond to U = 0, 2, and 4; data are for At = 0.125, 
except those with crossed symbols. Error bars result from 
uncertainties in the extrapolations ujm/2TTT — >■ (see text) 
and are only appreciable for p — 1. 



the smallest tti's (figure not shown); we then obtain the 
temperature-dependent Drude weight, D(T), appearing 
in Eq. (13). By varying the temperature, system size, 



and [/, we can generate plots of D(T), examples of which 
are shown in Figs. [8] and |9] As shown in Fig. [8J for 
a 12 X 12 lattice at half filling, in the non-interacting 
case the Drude weight clearly extrapolates to a non-zero 
value as T — > 0. For U > 0, D{T) vanishes at some 
temperature Ti^{L, U), which increases with U for a given 
L. Data for half filling in Fig. [9] show that at a fixed 
temperature, the Drude weight vanishes as the lattice size 
increases; that is, the points below Tq in Fig. [8] should 
approach the D = line for sufficiently large L. 

Away from half filling, the minus-sign problem pre- 
vents us from analyzing the size-dependence at very low 
temperatures, and we are restricted to data for /3 = 16 
for the densities p = 0.42, and 0.66, while keeping U < 2. 
Nonetheless, some important conclusions can be drawn 
from our analyses of the data for D{T) on finite-sized lat- 
tices: (1) we have found no evidence of a vanishing Drude 
weight at fixed, finite temperatures in the limit L — ^ oo, 
as previously suggested for the one-dimensional casepH 
(2) the dependence of D with 1/L, for fixed both tem- 
perature and on-site repulsion, is rather weak, without 
suffering from closed-shell effects, thus rendering extrap- 
olations towards L — > oo trustworthy. Once again, data 
for At = 0.125 are the same as those for At = 0.0625, 
within error bars. 

Our results therefore show that the Drude weight has 
been hitherto unjustifiably overlooked as a reliable probe 
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FIG. 10. (Color online) Dependence of the (normalized) 
Drude weight with the square of the 'time' interval, at fixed 
temperature and lattice size at half filling (triangles), and at 
p = 0.42 (circles). Error bars result from uncertainties in the 
extrapolations LUrn/^nT — >■ (see text). 



of the metal- insula! or transition; its use should be more 
widespread, given that it is free from closed-shell effects, 
and its clear-cut 'temperature-dependence' allows for an 
unambiguous characterization of insulating states. 



VI. SINGLE-PARTICLE EXCITATION GAP 

Another quantity used to infer the transport properties 
of the system is the single-particle excitation gap Asp(q), 
which is the minimum energy necessary to extract one 
fermion from the system, and is, essentially, related to 
the gap measurable in photoemission experiments. It can 
be obtained from the imaginary-time-dependent Green's 
function in reciprocal space, which for large r decays ex- 
ponentially, i.e., G{(Ip,t) ~ e~^=''(''^^'^ (see, e.g., Ref. 
[5^ . We can therefore obtain A^p through fits of QMC 
data for the Green's function, calculated at the Fermi 
wavevector for the electronic densities of interest. Figure 
|11| shows the imaginary-time dependence of the Green's 
function for the half-filled case. In the upper panel, the 
absence of a decay in the non-interacting case is a signa- 
ture of a metallic state, while the exponential decay in 
the lower panel results from a finite gap. The inset in (b) 
compares data obtained for two values of At: the time- 
dependent Green's functions lie on the same exponential 
curve, which illustrates that this quantity is also negligi- 
bly dependent on the At used. The size dependence of 
the gap is shown in Fig. [T2j for different values of U; for 
U = 2, one also sees that data for a smaller At lie on 
the same curve. The limiting (i.e., L — > oo) value of Asp 
increases from zero with increasing U, as expected; it is 
again clear that the value of At does not influence this 




FIG. 11. (Color online) Log-linear plot of the imaginary-time 
dependence of the Green's function G(q, r) at the Fermi wave 
vector qF, for different lattice sizes at half filling, p = 1.0, for 
the non-interacting (a), and interacting (b) cases. The error 
bars in (b) are due to statistical errors from averaging over 
different realizations, and equivalent qf points; here Ar = 
0.125. The inset includes data obtained with Ar — 0.0625, 
denoted by the corresponding crossed symbols from the main 
panel. 



extrapolation procedure. 

Figure [13] shows data for the Green's function for the 
density p — 0.42. In the non-interacting case, and dis- 
carding the data for L = 10, we see that the slope 
decreases as L increases, leading to a vanishing gap as 
i — > cxD, as one would expect for a metallic system; the 
data for i = 10 are completely off the mark, again as 
a result from the closed-shell density for this L. For 
the interacting case [Fig. 13 'b)], the Green's function for 



L 7^ 10 behaves in a way similar to that for the free case; 
again, the L = 10 case behaves completely differently 
from the others, bearing a negative gap as the signature 
of the closed-shell problem. In this respect, it is inter- 
esting to have in mind that the single-particle excitation 
gap provides a very clear indication that a closed-shell 
incident is at play for a given combination of p and L. 
For completeness, we note that, similarly to half filling 
(Fig. Ill, the dependence with At is negligible. 
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FIG. 12. (Color online) Finite-size dependence of the single 
particle excitation gap Asp(q^) at half- filling, for different 
values of the on-site repulsion. The error bars are due to the 
exponential fits to the data for G(qj^, r); see text. The crossed 
symbols denote the corresponding data for Ar=0.0625. 

VII. THE MINUS-SIGN PROBLEM 

In the present formulation of the QMC method, once 
the fermionic degrees of freedom are traced out, the role 
of Boltzmann factor in the partition function is played by 
the product of two determinants; see, e.g., Refs.lUIBl and 
[71 Since one cannot guarantee that this product is posi- 
tive definite for each configuration of the auxiliary fields, 
the averages are carried out in the ensemble of positive 
Boltzmann weights, at the expense of having to divide 
these averages by the average sign of the product of deter- 
minants, {sign). Therefore, when (sign) becomes signifi- 
cantly smaller than 1, the average values of most quanti- 
ties of interest become meaningless: this is the infamous 
'minus-sign problem'. It should be noted that other im- 
plementations of the QMC method also run into similar 
problems; see, e.g., Ref. 33|i 

This problem has eluded a variety of attempts of solu- 
tion proposed over the years; see, e.g., Ref. 7 for a partial 
list of references. For instance, once realized that simply 
ignoring the negative sign leads to serious discrepancies)^ 
attempts to use different Hubbard-Stratonovich trans- 
formations turned out to be fruitlessp^Ml the minus- 
sign problem has been alleviated with implementations 
of QMC constraining the sampling process, '^^'^^' from 
which a ground-state wave function is obtained. Other 
frameworks have been proposed to improve the sign 
problem,S2H321 fj^t systematic implementations compar- 
ing results for, e.g., correlation functions in the Hubbard 
model are, as far as we know, still unavailable. 

More recently, arguments have been giverP^l suggest- 
ing that there is no generic solution to the sign problem; 
instead, in the most favorable scenario, one may find spe- 
cial solutions for specific models.*^ 

In view of this, it is imperative to gather as much in- 
formation as possible about (sign). With this in mind. 




FIG. 13. (Color online) Same as Fig. |11[ but now for the 
electronic density p — 0A2. 

we define a quantity, k = 1 ~ p^n, directly related to 



the compressibility k defined in Sec. |III| Figure [14] shows 
that k reaches the value 1 at the densities correspond- 
ing to 'closed-sheir, as already discussed. In the same 
figure we also show {sign) as a function of p: interest- 
ingly, we see that it tracks k, in the sense that, at least 
for C/ < 2, it is harmless at densities such that k « 1 
(k w 0), but it can be seriously deleterious to the QMC 
averaging process when the system is more compressible, 
especially at larger values of U. Since a larger compress- 
ibility, in turn, corresponds to stronger density fluctua- 
tions, one may conclude that these are inherently linked 
with the minus-sign problem. It is worth noticing that 
improvements on co nvergence have been achieved within 
both projectoi'^SH g^j^^j fixed-node^ QMC simulations 
if closed-shell configurations are used as initial states; 
in addition, in Ref. HH it was also pointed out that the 
choice of closed-shell initial states led to larger {sign) 
than when open-shell initial states were taken. On the 
other hand, shell effects have also disrupted the density 
dependence need ed in the search for phase separation 
in the t-J modeL^^'^ Thus, while indications of an in- 
terplay between closed-shell and the minus-sign problem 
have been suggested in the past, Fig. [14] presents the first 
systematic evidence of this connection. 

It is also instructive to examine the behavior of {sign) 
with At. Figure 15 compares data for one lattice size. 




0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



FIG. 14. (Color online) Average sign (circles) of the fermionic determinant and R [(green) thick line for [/ = 0, and (green) 
triangles for f/ = 2; see text for the definition of k] as functions of electronic density, for different system sizes. Filled, half-filled, 
and empty circles respectively denote U — 2, 3 and 4; data for U — are with j3 = 30, while for [/ 7^ data are with /? — 16. 
For the sake of clarity, error bars were omitted, since they are smaller than data points. 



L = 10, but for different values oiU, (3, and the chemical 
potential /i. For /3 = 10, we see that for the closed shell 
configuration, p = 0.42, (sign) « 1 for all At in the 
range considered, for both U — 2 and 4; this feature 
is maintained when /3 is increased to 16, illustrating the 
harmlessness of (sign) at the closed-shell density. Doping 
slightly away, e.g., for an open shell configuration with 
p ~ 0.5, (sign) remains almost independent of Ar for 
U = 2, but acquires a significant dependence for U — A, 
leading to low values for small Ar; for (3 = 16, (sign) w 
for all At in the relevant range. Worse still, for [/ = 4 
and p ~ 0.7, (sign) is very close to zero for all values of 
At considered, for both f3 — 10 and 16; this should not 
come as a surprise, since Fig. [M] shows that k vanishes 
near this density for the 10 x 10 lattice. 



In Figure 16 we display the dependence of two aver- 
age local quantities with At^ , for two fixed values of the 
chemical potential, and for [/ = 4 and /? = 10. Notwith- 
standing the fact that systematic errors of order At^ are 



expected as a result of the Suzuki- Trotter decomposition, 
panel (a) shows that for p w —3.2, the proportionality 
constant is quite small, so that p — 0.42 over the whole 
range of At; in the open shell case, for which (sign) dete- 
riorates with decreasing At (see Fig. 15 1, the dependence 
of p with At^ is noticeable. By contrast, the lower panel 
shows that the double occupancy. 



{nnnii) 



(14) 



follows the expected linear dependence with At^ in both 
cases. This indicates that whenever {sign) is strongly de- 
pendent on At, one can still obtain meaningful averages 
by using solely the data for the largest values of At to 
extrapolate towards At = 0; though with less confidence, 
the same procedure could be adopted for /3 = 16. 
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FIG. 15. (Color online) Average sign of the fermionic de- 
terminant as a function of the square of Suzuki- Trotter 'time' 
interval, for a 10 x 10 lattice, with (a) /3 — 10, and (b) /3 — 16. 
Black squares and (green) up triangles respectively corre- 
spond to the closed-shell density p = 0.42, and to p « 0.5; 
half-filled and filled symbols respectively correspond to [/ = 2 
and 4. 
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FIG. 16. (Color online) Dependence of average values of 
(a) electronic density, and (b) double occupancy (see text) 
with the square of the Suzuki-Trotter 'time' interval, for two 
values of the chemical potential; U, /3, and lattice size are 
fixed. The extrapolated values, obtained from the fitting of a 
straight line through all points, are shown in red at At — 0. 



ulations at finite temperatures for the homogeneous Hub- 
bard model on the square lattice, and commonly used to 
locate insulating behavior. Our results show that 'closed- 
shell' effects, which introduce important (though artifi- 
cial) gaps in the spectrum, may lead to false insulating 
behavior of the compressibility, of the conductivity, and 
of the charge gap at certain combinations of occupation 
and linear lattice size, L; in situations in which a long 
series of lattice sizes cannot be obtained, this may jeop- 
ardize extrapolations towards L — )• c». We have also 
assessed corrections to the dc-conductivity, which are ne- 
glected when a Laplace transform is avoided through a 
simplifying prescription, and found that the latter is not 
generically valid due to the absence of a sufficiently small 
energy scale in the problem; though quite appealing, fit- 
tings to experimental data with the conductivity thus 
obtained should be avoided. The Drude weight, on the 
other hand, suffers from more controllable finite-size and 
finite-temperature effects. At half filling, and at a fixed 
low temperature, it vanishes with a power law in 1/L, 
the exponent of which depends on U; away from half 
filling, the Drude weight is only weakly dependent on ei- 
ther temperature and system size, being free from the 
spurious behavior found in other quantities. Therefore, 
amongst all quantities discussed here, the Drude weight 
is certainly the most reliable one to use in situations for 
which the data are limited to a restricted set of system 
sizes. 

In addition, we have also presented numerical evidence 
showing that the sign of the fermionic determinant tracks 
the compressibility: for densities at which the system is 
'incompressible', as a result of a gap due to the finiteness 
of the lattice, (sign) w 1, at least for U < 2. However, in- 
between two successive 'incompressible' densities, (sign) 
deteriorates steadily as U increases. This behavior is 
suggestive that strong density fluctuations may be linked 
to the 'minus-sign problem'. We have also investigated 
the influence of the imaginary-time interval At on the 
behavior of {sign) and of some (local) average quantities. 

All analysed quantities can be fitted to a linear de- 
pendence with At^ , as expected from the Suzuki- Trotter 
discretization, although at the closed-shell density, the 
slopes for both (sign) and the density p are very small. 
The At^ dependence is indicative that for some densi- 
ties, one can confidently use data for 'large' At (i.e., 
those leading to {sign) ^ 0.5) to perform extrapolations 
(towards At — >■ 0) of average values. 



VIII. CONCLUSION 

In conclusion, we have thoroughly examined the be- 
havior of several quantities, obtained through QMC sim- 
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